package edu.hawaii.senin.iclouds.simple;

import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Calendar;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Observable;
import java.util.Observer;
import java.util.TreeMap;

import org.apache.commons.math.random.RandomDataImpl;
import org.jdom.Document;
import org.jdom.Element;
import org.jdom.JDOMException;
import org.jdom.input.SAXBuilder;

import edu.hawaii.senin.iclouds.basics.Atom;
import edu.hawaii.senin.iclouds.basics.Constants;
import edu.hawaii.senin.iclouds.basics.Interval;
import edu.hawaii.senin.iclouds.basics.Molecule;
import edu.hawaii.senin.iclouds.basics.Reaction;
import edu.hawaii.senin.iclouds.basics.Specie;
import edu.hawaii.senin.iclouds.controller.Controller;
import edu.hawaii.senin.iclouds.util.CloudPhysicalParameters;
import edu.hawaii.senin.iclouds.util.CloudPopulation;
import edu.hawaii.senin.iclouds.util.CloudReactions;

/**
 * This class represents simple model of interstellar cloud.
 * 
 * @author Pavel V. Senin
 * 
 */
public class Cloud extends Observable implements Runnable, Observer {

  /**
   * Grain size.
   */
  private Double cloudGrainSize = 1000.0;

  /**
   * Contains all atoms and molecules within the cloud along with abundances.
   */
  private HashMap<String, Specie> cloudPopulation;
  /**
   * Contains all atoms and molecules within the grain surface along with their abundances.
   */
  private HashMap<String, Specie> grainPopulation;
  /**
   * Contains all specified reactions for that gaseous phase along with rates.
   */
  private ArrayList<Reaction> reactions;

  /**
   * Specifies temperature within this cloud.
   */
  private Double cloudTemperature;

  /**
   * Specifies density within this cloud.
   */
  private Double cloudDensity;

  /**
   * Specifies current time frame.
   */
  private Double time;

  /**
   * Specifies accretion time.
   */
  private Double cloudAccretionTimeSpan;

  /**
   * specifies initial cloud configuration.
   */
  private Document configFileName;

  /**
   * Limit of iterations, will exit from main loop if even not converge yet.
   */
  private Integer iterationsLimit;

  /**
   * Number of reaction sites?
   */
  private Integer cloudGrainsNumber;

  /**
   * Barrier width.
   */
  public Double width;

  /**
   * Counts iterations.
   */
  private Integer iterationCounter;

  /**
   * Seed for random generator.
   */
  private Integer randomSeed;

  /**
   * Random number generator handler.
   */
  private RandomDataImpl randomGenerator;

  /**
   * How many species accreted at once.
   */
  private Integer cloudAccretionBatchSize;

  private Calendar calendar = Calendar.getInstance();

  private boolean ready;

  private BufferedWriter fileWriter;

  private boolean stopSimulation;

  /**
   * Instantiates and initializes cloud model. Calculates all needed data from supplied by
   * configuration initial data. Does exactly next moves:
   * <ul>
   * <li>fixes local time and iteration;</li>
   * <li>sets velocities;</li>
   * <li>sets activation barrier width;</li>
   * <li>fixes all reactions accretion rate;</li>
   * </ul>
   * 
   */
  public void initialize() {

    // init random generator with the same seed.
    this.randomSeed = 112;
    this.randomGenerator = new RandomDataImpl();
    this.randomGenerator.reSeed(this.randomSeed);
    this.grainPopulation = new HashMap<String, Specie>();

    // check and set the default time
    if (null == this.time) {
      this.setTime(0.0);
    }

    // check and set this iteration counter
    if (null == this.iterationCounter) {
      this.setIterationCtr(0);
    }

    // check and set velocities for species
    for (Specie sp : this.cloudPopulation.values()) {
      sp.setVelocity(Math.sqrt(3 * Constants.boltzmann * this.cloudTemperature
          / (Constants.atomicMassUnit * sp.getAmu())));
    }

    // set-up the barrier width
    this.width = 1E-10;

    // Calculate the probability of overcoming an energy barrier
    for (Reaction r : this.reactions) {

      String type = r.getType();

      /*---------------------------------------------*/
      if (type.equalsIgnoreCase("QT")) {
        // (type[j] EQ 'QT') : chi[j] = nu*exp(-((2*a/h_bar)*(sqrt(2.*mass_H*k*E[j])))
        r.setChi(Constants.verticalOscillation
            * Math.exp((-2 * this.width / Constants.reducedPlank)
                * Math.sqrt(2 * Constants.hydrodgeneMass * Constants.boltzmann * r.getEBarrier())));
      }
      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF")) && (r.getReactants().contains("H"))) {
        // (type[j] EQ 'DF') AND (Reactant[j] EQ 'H'): chi[j] =
        // nu*exp(-((2*a/h_bar)*(sqrt(2.*mass_H*k*E[j]))))
        r.setChi(Constants.verticalOscillation
            * Math.exp(-((2 * this.width / Constants.reducedPlank) * Math.sqrt(2
                * Constants.hydrodgeneMass * Constants.boltzmann * r.getEBarrier()))));
      }
      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF")) && (r.getReactants().contains("D"))) {
        // (type[j] EQ 'DF') AND (Reactant[j] EQ 'H'): chi[j] =
        // nu*exp(-((2*a/h_bar)*(sqrt(2.*mass_H*k*E[j]))))
        r.setChi(Constants.verticalOscillation
            * Math.exp(-(r.getEBarrier() / this.cloudTemperature)));
      }
      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF")) && (r.getReactants().contains("O"))) {
        // (type[j] EQ 'DF') AND (Reactant[j] EQ 'O'): chi[j] = nu*exp(-(240/T))
        r.setChi(Constants.verticalOscillation * Math.exp(-(240 / this.cloudTemperature)));
      }

      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF_O")) && (r.getReactants().contains("O"))) {
        // (type[j] EQ 'DF_O') AND (Reactant[j] EQ 'O'): chi[j] = nu*exp(-(E[j]/T))
        r.setChi(Constants.verticalOscillation
            * Math.exp(-(r.getEBarrier() / this.cloudTemperature)));
      }

      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF_C")) && (r.getReactants().contains("CO"))) {
        // (type[j] EQ 'DF_C') AND (Reactant[j] EQ 'CO'): chi[j] = nu*exp(-(E[j]/T))
        r.setChi(Constants.verticalOscillation
            * Math.exp(-(r.getEBarrier() / this.cloudTemperature)));
      }

      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF_N")) && (r.getReactants().contains("N"))) {
        // (type[j] EQ 'DF_N') AND (Reactant[j] EQ 'N'): chi[j] = nu*exp(-(E[j]/T))
        r.setChi(Constants.verticalOscillation
            * Math.exp(-(r.getEBarrier() / this.cloudTemperature)));
      }

      /*---------------------------------------------*/
      if ((type.equalsIgnoreCase("DF_C")) && (r.getReactants().contains("C"))) {
        // (type[j] EQ 'DF_C') AND (Reactant[j] EQ 'C'): chi[j] = nu*exp(-(E[j]/T))
        r.setChi(Constants.verticalOscillation
            * Math.exp(-(r.getEBarrier() / this.cloudTemperature)));
      }
    }

    // // calculate accretion time NOT USED SO FAR
    // this.cloudAccretionTimeSpan = 0.0;
    // for (Specie sp : this.cloudPopulation.values()) {
    // this.cloudAccretionTimeSpan += sp.getRate(this.cloudGrainSize);
    // }
    // this.cloudAccretionTimeSpan = 1 / this.cloudAccretionTimeSpan;

  }

  private void setIterationCtr(int counterValue) {
    this.iterationCounter = counterValue;

  }

  /**
   * Sets cloud time.
   * 
   * @param time new time to set.
   */
  private void setTime(Double time) {
    this.time = time;

  }

  /**
   * Returns current cloud population.
   * 
   * @return actual population.
   */
  public ArrayList<Object> getPopulation() {
    ArrayList<Object> currentPopulation = new ArrayList<Object>();
    for (String name : this.cloudPopulation.keySet()) {
      currentPopulation.add(this.cloudPopulation.get(name));
    }
    return currentPopulation;
  }

  /**
   * Parses population section of the config.
   * 
   * @param species JDOM formatted list of all species.
   */
  private void parsePopulation(List species) {
    Iterator itr = species.iterator();
    while (itr.hasNext()) {

      Element el = (Element) itr.next();

      String name = el.getChild("name").getValue();
      String type = el.getChild("type").getValue();
      Double amu = Double.valueOf(el.getChild("amu").getValue());
      Double abundance = Double.valueOf(el.getChild("abundance").getValue());
      Double fraction = Double.valueOf(el.getChild("fraction").getValue());

      if ("atom".equalsIgnoreCase(el.getName())) {
        this.cloudPopulation.put(name, new Atom(abundance, -1.0, amu, fraction, name, type, -1.0,
            0L));
      }
      if ("molecule".equalsIgnoreCase(el.getName())) {
        this.cloudPopulation.put(name, new Molecule(abundance, -1.0, amu, fraction, name, type,
            -1.0, 0L));
      }
    }
  }

  /**
   * Parses reactions section of the config.
   * 
   * @param react JDOM formatted list of all reactions.
   */
  private void parseReactions(List react) {
    // <reaction>
    // <name>foo</name>
    // <rate>1.0</rate>
    // <reactants>
    // <atom>H</atom>
    // <atom>C</atom>
    // </reactants>
    // <products>
    // <molecule>CH</molecule>
    // </products>
    // <energy_barrier>105.0</energy_barrier>
    // <type>DF</type>
    // </reaction>
    Iterator itr = react.iterator();
    while (itr.hasNext()) {

      Element el = (Element) itr.next();

      String name = el.getChild("name").getValue();
      String type = el.getChild("type").getValue();
      Double barrier = Double.valueOf(el.getChild("energy_barrier").getValue());
      // Double rate = Double.valueOf(el.getChild("rate").getValue());

      ArrayList<String> reactants = getReactants(el.getChild("reactants").getChildren());
      ArrayList<String> products = getProducts(el.getChild("products").getChildren());

      // String nme = "";
      // for(String s : reactants){
      // nme = nme.concat(s);
      // }
      // nme = nme.concat("_");
      // for(String s : products){
      // nme = nme.concat(s);
      // }
      //      
      // System.out.println("<reaction>");
      // System.out.println("<name>" + nme + "</name>");
      // System.out.println("<rate>"+rate+"</rate>");
      // System.out.println("<reactants>");
      // for(String s : reactants){
      // if(s.length()==1){
      // System.out.println("<atom>"+s+"</atom>");
      // }else{
      // System.out.println("<molecule>"+s+"</molecule>");
      // }
      // }
      // System.out.println("</reactants>");
      // System.out.println("<products>");
      // for(String s : products){
      // if(s.length()==1){
      // System.out.println("<atom>"+s+"</atom>");
      // }else{
      // System.out.println("<molecule>"+s+"</molecule>");
      // }
      // }
      //      
      // System.out.println("</products>");
      // System.out.println("<energy_barrier>" + barrier + "</energy_barrier>");
      // System.out.println("<type>" + type + "</type>");
      // System.out.println("</reaction>");

      this.reactions.add(new Reaction(name, barrier, 1.0, type, reactants, products));
    }
  }

  /**
   * Extracts reactants from the list.
   * 
   * @param reactants list to extract from.
   * @return the text-element list.
   */
  private ArrayList<String> getReactants(List reactants) {

    ArrayList<String> rList = new ArrayList<String>();
    Iterator itr = reactants.iterator();
    while (itr.hasNext()) {
      rList.add(((Element) itr.next()).getValue());
    }

    return rList;
  }

  /**
   * Extracts reaction products from the list.
   * 
   * @param products list of products to extract from.
   * @return the text-element list.
   */
  private ArrayList<String> getProducts(List products) {

    ArrayList<String> pList = new ArrayList<String>();
    Iterator itr = products.iterator();
    while (itr.hasNext()) {
      pList.add(((Element) itr.next()).getValue());
    }

    return pList;
  }

  /**
   * Sets the clousd temperature.
   * 
   * @param temperature new temperature.
   */
  public void setTemperature(Double temperature) {
    this.cloudTemperature = temperature;

  }

  /**
   * Returns cloud temperature.
   * 
   * @return cloud temperature.
   */
  public Double getTemperature() {
    return this.cloudTemperature;
  }

  /**
   * Sets the cloud density.
   * 
   * @param newDensity new density.
   */
  public void setDensity(Double newDensity) {
    this.cloudDensity = newDensity;

  }

  /**
   * Returns cloud density.
   * 
   * @return cloud density.
   */
  public Double getDensity() {
    return this.cloudDensity;
  }

  /**
   * Returns current cloud snapshot in Plain text format.
   * 
   * @param fileName file name to dump to.
   * 
   */
  public void writeSnapshot() {

    if (null == this.fileWriter) {
      assert true;
    }
    else {
      try {
        this.fileWriter.write("Iterations left: " + String.valueOf(this.iterationsLimit));
        this.fileWriter.write("\nTime:        " + String.valueOf(this.getTime()));
        this.fileWriter.write("\n\nN sites:     " + String.valueOf(this.cloudGrainsNumber));
        this.fileWriter.write("\nDensity:     " + String.valueOf(this.getDensity()));
        this.fileWriter.write("\nTemperature: " + String.valueOf(this.getTemperature()));
        this.fileWriter.write("\ngrainSize:   " + String.valueOf(this.cloudGrainSize));
        // this.fileWriter.write("\nActivation H2CO: " + String.valueOf(this.activation_H2CO));
        // this.fileWriter.write("\nactivation_H2CO_computed: " +
        // String.valueOf(this.activation_H2CO_computed));
        this.fileWriter.write("\naccretion time computed : "
            + String.valueOf(this.cloudAccretionTimeSpan));
        this.fileWriter.write("\n\nCurrent Population: \n");
        for (Specie sp : this.cloudPopulation.values()) {
          this.fileWriter.write(" atom: " + sp.getName() + ", amu: " + String.valueOf(sp.getAmu())
              + ", abun: " + String.valueOf(sp.getAbundance()) + ", fraction: "
              + String.valueOf(sp.getFraction()) + ", vel: " + String.valueOf(sp.getVelocity())
              + ", acc.rate: " + String.valueOf(sp.getRate(this.cloudGrainSize)) + "\n");
        }
        this.fileWriter.write("\nReactions configured: \n");
        for (Reaction r : this.reactions) {
          this.fileWriter.write("  Reaction: " + r.getName() + ", type: " + r.getType() + ", Chi: "
              + String.valueOf(r.getChi()) + ", rate: " + String.valueOf(r.getRate())
              + ", Ebarrier: " + String.valueOf(r.getEBarrier()) + ", reactants: "
              + r.getReactants().toString() + ", products: " + r.getProducts().toString() + "\n");
        }
        this.fileWriter.write("\nAccretion probabilities: \n");
        Double sum = 0.0;
        for (Specie sp : this.cloudPopulation.values()) {
          this.fileWriter.write("  atom: " + sp.getName() + ", accretion probability: "
              + String.valueOf(sp.getAccretionProbability()) + "\n");
          sum = sum + sp.getAccretionProbability();

        }
        this.fileWriter.write("Total (sum) is " + String.valueOf(sum) + "\n");
        // this.fileWriter.write(" For CO_OH probability is: " + this.chiCOOH + "\n");
        // this.fileWriter.write(" For H2OH probability is: " + this.chiH2OH + "\n");
        this.fileWriter.flush();
      }
      catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }
    }

  }

  /**
   * Gets current cloud time.
   * 
   * @return time in days.
   */
  private Double getTime() {
    return this.time;
  }

  /**
   * Updates all accretion probabilities.
   */
  public void updateAccretionRates() {
    // Original code looks like:
    // n_abundance = (10.^iterate)*1L
    // n_O2 = (f_O2*(A_O-(ABS((1-f_C)*A_C)))*n_abundance)*mo2 ;abundance of molecular oxygen
    // n_O = (1-2*f_O2)*(ABS(A_O-((1-f_C)*A_C)))*n_abundance ;abundance of oxygen
    // n_N = (A_N*f_N)*n_abundance
    // n_N2 = (((1-f_N)*A_N)/2)*n_abundance
    // n_C = (A_C*f_C)*n_abundance
    // n_CO = ((1-f_C)*A_C)*n_abundance ;abundance of carbon monoxide

    for (String specie : this.cloudPopulation.keySet()) {
      if (specie.equals("O")) {
        // n_O = (1-2*f_O2)*(ABS(A_O-((1-f_C)*A_C)))*n_abundance ;abundance of oxygen
        Atom carbon = (Atom) this.cloudPopulation.get("C");
        Atom oxygene = (Atom) this.cloudPopulation.get("O");
        Molecule oxygeneM = (Molecule) this.cloudPopulation.get("O2");
        Double ap = (1 - oxygeneM.getFraction() * 2)
            * Math.abs(oxygene.getAbundance() - (1 - carbon.getFraction()) * carbon.getAbundance())
            * this.cloudDensity;
        oxygene.setAccretionProbability(ap);
      }
      if (specie.equals("O2")) {
        // n_O2 = (f_O2*(A_O-(ABS((1-f_C)*A_C)))*n_abundance)*mo2 ;abundance of molecular oxygen
        Atom carbon = (Atom) this.cloudPopulation.get("C");
        Atom oxygene = (Atom) this.cloudPopulation.get("O");
        Double ap = (oxygene.getFraction() * (oxygene.getAbundance() - Math.abs(1 - carbon
            .getFraction())
            * carbon.getAbundance()))
            * this.cloudDensity;
        ((Molecule) this.cloudPopulation.get("O2")).setAccretionProbability(ap);
      }
      if (specie.equals("N")) {
        // n_N = (A_N*f_N)*n_abundance
        Atom nitrogene = (Atom) this.cloudPopulation.get("N");
        Double ap = (nitrogene.getAbundance() * nitrogene.getFraction()) * this.cloudDensity;
        nitrogene.setAccretionProbability(ap);
      }
      if (specie.equals("N2")) {
        // n_N2 = (((1-f_N)*A_N)/2)*n_abundance
        Molecule nitrogeneM = (Molecule) this.cloudPopulation.get("N2");
        Atom nitrogene = (Atom) this.cloudPopulation.get("N");
        Double ap = ((1 - nitrogene.getFraction()) * nitrogene.getAbundance() / 2)
            * this.cloudDensity;
        nitrogeneM.setAccretionProbability(ap);
      }
      if (specie.equals("C")) {
        // n_C = (A_C*f_C)*n_abundance
        Atom carbon = (Atom) this.cloudPopulation.get("C");
        Double ap = (carbon.getAbundance() * carbon.getFraction()) * this.cloudDensity;
        carbon.setAccretionProbability(ap);
      }
      if (specie.equals("CO")) {
        // n_CO = ((1-f_C)*A_C)*n_abundance ;abundance of carbon monoxide
        Atom carbon = (Atom) this.cloudPopulation.get("C");
        Molecule carbonOxide = (Molecule) this.cloudPopulation.get("CO");
        Double ap = (1 - carbon.getFraction()) * carbon.getAbundance() * this.cloudDensity;
        carbonOxide.setAccretionProbability(ap);
      }
      if (specie.equals("H")) {
        // this is a constant
        Atom hydrogen = (Atom) this.cloudPopulation.get("H");
        hydrogen.setAccretionProbability(1.0);
      }
      if (specie.equals("D")) {
        // this is a constant
        Atom deiterium = (Atom) this.cloudPopulation.get("D");
        Double ap = deiterium.getFraction() * deiterium.getAbundance() * this.cloudDensity;
        deiterium.setAccretionProbability(ap);
      }
    }

    Double normalizer = 0.0;
    for (Specie sp : this.cloudPopulation.values()) {
      normalizer = normalizer + sp.getAccretionProbability() * sp.getVelocity();
    }
    for (Specie sp : this.cloudPopulation.values()) {
      sp.setAccretionProbability((sp.getAccretionProbability() * sp.getVelocity()) / normalizer);
    }
  }

  /**
   * Reports number of sites on grain surface.
   * 
   * @return number of sites.
   */
  public Integer getNSites() {
    return this.cloudGrainsNumber;
  }

  /**
   * Reports cloud density.
   * 
   * @return cloud density.
   */
  public Double getAbundance() {
    return this.cloudDensity;
  }

  /**
   * Reports limit of accreted at once species.
   * 
   * @return limit of accreted species.
   */
  public Integer getAccretionStep() {
    return this.cloudAccretionBatchSize;
  }

  /**
   * Reports grain current grain size.
   * 
   * @return grain size.
   */
  public Double getGrainSize() {
    return this.cloudGrainSize;
  }

  /**
   * Reports defined iteration limit.
   * 
   * @return iteration limit.
   */
  public Integer getIterationLimit() {
    return this.iterationsLimit;
  }

  /**
   * Does accretion step from gas phase to grain surface.
   * 
   * @return map where each key is the name of specie and value is number of accreted species.
   */
  public TreeMap<String, Specie> cycleAccretion() {
    /**
     * Map that holds all species to accrete. The interval is mapped to certain specie.
     */
    TreeMap<String, Interval> species = new TreeMap<String, Interval>();
    /**
     * Map that "accumulates" species. Maps accreted volume to specie name.
     */
    TreeMap<String, Specie> accretionSample = new TreeMap<String, Specie>();

    // initialize all intervals according the rates
    Double currentLow = 0.0;
    for (Specie sp : this.cloudPopulation.values()) {
      if (sp.getAccretionProbability().equals(0D)) {
        // System.out.println("omitted" + sp.getName());
      }
      else {
        species.put(sp.getName(), new Interval(currentLow, currentLow
            + sp.getAccretionProbability()));
        currentLow = currentLow + sp.getAccretionProbability();
        accretionSample.put(sp.getName(), sp.clone());
      }
    }
    // accrete all guys
    for (int i = 0; i < this.cloudAccretionBatchSize; i++) {
      Double rnd = this.randomGenerator.nextUniform(0D, 1D);
      for (String sp : species.keySet()) {
        if (species.get(sp).within(rnd)) {
          accretionSample.get(sp).incrementAccumulatorValue(1L);
          break;
        }
      }
    }

    return accretionSample;
  }

  private String getTimeSpan(Long startTime, Long endTime) {
    Long delta = endTime - startTime;
    Long hours = delta / 3600000;
    Long minutes = (delta % 3600000) / 60000;
    Long seconds = (delta - hours * 3600000 - minutes * 60000) / 1000;
    Long milliseconds = delta - hours * 3600000 - minutes * 60000 - seconds * 1000;
    return hours + "h " + minutes + "m " + seconds + "s " + milliseconds + " ms";
  }

  /**
   * Gets over reaction network with provided accreted species.
   * 
   * @param newComers newly accreted species.
   * @return new grain surface picture.
   */
  public void SurfaceCycle(TreeMap<String, Specie> newComers) {

    boolean inLoop = true;
    // cleanUp newComers
    while (inLoop) {
      inLoop = false;
      for (Specie sp : newComers.values()) {
        if (sp.getAccumulatorValue().equals(0L)) {
          newComers.remove(sp.getName());
          inLoop = true;
          break;
        }
      }
    }

    // main cycle
    inLoop = true;

    // Map that holds all species to accrete. The interval is mapped to certain specie.
    TreeMap<String, Interval> chosenIntervals = new TreeMap<String, Interval>();
    // holds all reactions for chosen specie
    TreeMap<String, Reaction> possibleReactions = new TreeMap<String, Reaction>();
    TreeMap<String, Reaction> coreactants = new TreeMap<String, Reaction>();

    Integer iteration = 0;

    Long intervalsStartTime, intervalsEndTime, intervalsTotalTime = 0L;
    Long possibleReactionsStartTime, possibleReactionsEndTime, possibleReactionsTotalTime = 0L;

    while (inLoop) {

      possibleReactions.clear();
      coreactants.clear();

      // choose the victim!
      intervalsStartTime = System.currentTimeMillis();
      chosenIntervals = refreshIntervals(chosenIntervals, newComers);
      Double rnd = this.randomGenerator.nextUniform(0D, 1D);
      String choosenBeing = chooseBeing(rnd, chosenIntervals);
      intervalsEndTime = System.currentTimeMillis();
      intervalsTotalTime += (intervalsEndTime - intervalsStartTime);

      // go over it, we are having all possible reactions after this cycle
      // basically we are filtering all reactions with choosenBeing at first and all rections that
      // possible due to the presence of co-reactants at second

      if (this.cloudPopulation.get(choosenBeing).isActive()) {

        //active - choosing 
        possibleReactionsStartTime = System.currentTimeMillis();
        possibleReactions = getPossibleReactions(choosenBeing, possibleReactions);
        possibleReactionsEndTime = System.currentTimeMillis();
        possibleReactionsTotalTime += (possibleReactionsEndTime - possibleReactionsStartTime);

        if (possibleReactions.size() == 0) {
          if (choosenBeing.equalsIgnoreCase("H")) {
            newComers.get("H").decrementAccumulatorValue(1L);
            if (newComers.get("H").getAccumulatorValue() == 0) {
              newComers.remove("H");
            }
          }
          else {
            if (this.grainPopulation.containsKey(choosenBeing)) {
              this.grainPopulation.get(choosenBeing).incrementAccumulatorValue(1L);
            }
            else {
              Specie clone = newComers.get(choosenBeing).clone();
              clone.setAccumulatorValue(1L);
              this.grainPopulation.put(clone.getName(), clone);
            }
            if (newComers.get(choosenBeing).getAccumulatorValue() == 1) {
              newComers.remove(choosenBeing);
            }
            else {
              newComers.get(choosenBeing).decrementAccumulatorValue(1L);
            }
          }// if not H

        }
        else {
          // possible reactions are not zero - go for chemistry there
          if (choosenBeing.equalsIgnoreCase("H")) {
            for (String coreactant : possibleReactions.keySet()) {
              coreactants.put(coreactant, possibleReactions.get(coreactant));
            }
            TreeMap<String, Interval> intervals = getHIntervalsSurface(coreactants);

            // find coreactant
            rnd = this.randomGenerator.nextUniform(0D, 1D);
            Reaction choosenReaction = null;
            String coreactant = null;
            for (String sp : intervals.keySet()) {
              if (intervals.get(sp).within(rnd)) {
                choosenReaction = possibleReactions.get(sp);
                coreactant = sp;
                break;
              }
            }
            // do H accounting
            if (newComers.get("H").getAccumulatorValue() == 1) {
              newComers.remove("H");
            }
            else {
              newComers.get("H").decrementAccumulatorValue(1L);
            }
            // remove coreactant from grain
            if (this.grainPopulation.get(coreactant).getAccumulatorValue() == 1) {
              this.grainPopulation.remove(coreactant);
            }
            else {
              this.grainPopulation.get(coreactant).decrementAccumulatorValue(1L);
            }
            // add product to surface
            ArrayList<String> productSpecies = choosenReaction.getProducts();
            for (String productName : productSpecies) {
              if (this.grainPopulation.containsKey(productName)) {
                this.grainPopulation.get(productName).incrementAccumulatorValue(1L);
              }
              else {
                this.grainPopulation.put(productName, new Molecule(0D, 0D, 0D, 0D, productName,
                    "passive", 0D, 1L));
              }
            }// for
          }
          else {
            // here goes something that not H - different algorithm for interval
            // collect all possible coreactants and their abundances.
            for (String coreactant : possibleReactions.keySet()) {
              coreactants.put(coreactant, possibleReactions.get(coreactant));
            }
            TreeMap<String, Interval> intervals = getIntervalsSurface(coreactants);

            // find coreactant
            rnd = this.randomGenerator.nextUniform(0D, 1D);
            Reaction choosenReaction = null;
            String coreactant = null;
            for (String sp : intervals.keySet()) {
              if (intervals.get(sp).within(rnd)) {
                choosenReaction = possibleReactions.get(sp);
                coreactant = sp;
                break;
              }
            }
            // do accounting
            if (newComers.get(choosenBeing).getAccumulatorValue() == 1) {
              newComers.remove(choosenBeing);
            }
            else {
              newComers.get(choosenBeing).decrementAccumulatorValue(1L);
            }
            // remove coreactant from grain
            if (this.grainPopulation.get(coreactant).getAccumulatorValue() == 1) {
              this.grainPopulation.remove(coreactant);
            }
            else {
              this.grainPopulation.get(coreactant).decrementAccumulatorValue(1L);
            }
            // add product to surface
            ArrayList<String> productSpecies = choosenReaction.getProducts();
            for (String productName : productSpecies) {
              if (this.grainPopulation.containsKey(productName)) {
                this.grainPopulation.get(productName).incrementAccumulatorValue(1L);
              }
              else {
                this.grainPopulation.put(productName, new Molecule(0D, 0D, 0D, 0D, productName,
                    "passive", 0D, 1L));
              }
            }// for
          }
        }// possible reactions are not empty
      }// specie is active
      else {
        // not active - just accrete and collect
        if (this.grainPopulation.containsKey(choosenBeing)) {
          this.grainPopulation.get(choosenBeing).incrementAccumulatorValue(1L);
        }
        else {
          Specie clone = newComers.get(choosenBeing).clone();
          clone.setAccumulatorValue(1L);
          this.grainPopulation.put(clone.getName(), clone);
        }
        if (newComers.get(choosenBeing).getAccumulatorValue() == 1) {
          newComers.remove(choosenBeing);
        }
        else {
          newComers.get(choosenBeing).decrementAccumulatorValue(1L);
        }
      }
      if (newComers.size() == 0) {
        inLoop = false;
      }
      else if (newComers.size() == 1) {
        if (null != newComers.get("H")) {
          inLoop = false;
        }
      }
      iteration++;
    }

  }

  private String chooseBeing(Double rnd, TreeMap<String, Interval> chosenIntervals) {
    String res = null;
    for (String sp : chosenIntervals.keySet()) {
      // System.out.println(sp + " " + chosenIntervals.get(sp).getLowBound() + "-"
      // + chosenIntervals.get(sp).getHighBound());
      if (chosenIntervals.get(sp).within(rnd)) {
        res = sp;
        // break;
      }
    }
    return res;
  }

  private TreeMap<String, Reaction> getPossibleReactions(String choosenBeing,
      TreeMap<String, Reaction> possibleReactions) {
    possibleReactions.clear();
    for (Reaction r : this.reactions) {
      if (r.getReactants().contains(choosenBeing)) {
        ArrayList<String> cr = r.getReactants();
        for (String reactant : cr) {
          if (!reactant.equalsIgnoreCase(choosenBeing)) {
            if (this.grainPopulation.containsKey(reactant)) {
              possibleReactions.put(reactant, r);
            }
          }
        }
      }
    }
    return possibleReactions;
  }

  private TreeMap<String, Interval> getIntervalsSurface(TreeMap<String, Reaction> coreactants) {
    TreeMap<String, Interval> res = new TreeMap<String, Interval>();
    Double sum = 0D;
    // System.out.println(coreactants.keySet());
    for (String s : coreactants.keySet()) {
      // System.out.println(s);
      // System.out.println(this.grainPopulation.get(s).getAccumulatorValue().doubleValue());
      // System.out.print(coreactants.get(s).getName());
      // System.out.print(coreactants.get(s).getReactants());
      // System.out.print(coreactants.get(s).getProducts());
      // System.out.print(coreactants.get(s).getType());
      // System.out.println(coreactants.get(s).getChi());
      sum += this.grainPopulation.get(s).getAccumulatorValue().doubleValue()
          * coreactants.get(s).getChi();
    }
    Double low = 0D;
    for (String s : coreactants.keySet()) {
      Double high = low
          + (this.grainPopulation.get(s).getAccumulatorValue().doubleValue() * coreactants.get(s)
              .getChi()) / sum.doubleValue();
      res.put(s, new Interval(low, high));
      low = high;
    }
    return res;
  }

  private TreeMap<String, Interval> getHIntervalsSurface(TreeMap<String, Reaction> coreactants) {
    TreeMap<String, Interval> res = new TreeMap<String, Interval>();
    Long sum = 0L;
    for (String s : coreactants.keySet()) {
      sum += this.grainPopulation.get(s).getAccumulatorValue();
    }
    Double low = 0D;
    for (String s : coreactants.keySet()) {
      Double high = low + this.grainPopulation.get(s).getAccumulatorValue().doubleValue()
          / sum.doubleValue();
      res.put(s, new Interval(low, high));
      low = high;
    }
    return res;
  }

  /**
   * Updates intervals for a random choice of reacting species.
   * 
   * @param species space where put the map.
   * @return all intervals according current surface population.
   */
  private TreeMap<String, Interval> refreshIntervals(TreeMap<String, Interval> species,
      TreeMap<String, Specie> iPopulation) {
    species.clear();
    Long totalPopulation = 0L;
    for (Specie sp : iPopulation.values()) {
      totalPopulation += sp.getAccumulatorValue();
    }
    Double currentLow = 0.0;
    for (Specie sp : iPopulation.values()) {
      Double highBound = currentLow + sp.getAccumulatorValue().doubleValue()
          / totalPopulation.doubleValue();
      species.put(sp.getName(), new Interval(currentLow, highBound));
      currentLow = highBound;
    }
    return species;
  }

  /**
   * Reports current grain population.
   * 
   * @return grain population.
   */
  public HashMap<String, Specie> getGrainPopulation() {
    return this.grainPopulation;
  }

  /**
   * Loads xml formatted initial cloud configuration.
   * 
   * @param configName the filename of config.
   * @throws ReadConfigException if config file is not correctly formatted.
   */
  public void loadConfig(String configName) throws ReadConfigException {
    SAXBuilder builder = new SAXBuilder();
    this.configFileName = null;
    try {
      BufferedReader reader = new BufferedReader(new FileReader(new File(configName)));
      this.configFileName = builder.build(reader);
    }
    catch (FileNotFoundException e) {
      setChanged();
      notifyObservers("model : can't really open cloud config file" + configName + "\n"
          + e.getStackTrace());
    }
    catch (JDOMException e) {
      setChanged();
      notifyObservers("model : cloud config file is not well formed " + configName + "\n"
          + e.getStackTrace());
    }
    catch (IOException e) {
      setChanged();
      notifyObservers("model : IO exception while parsing config file  " + configName + "\n"
          + e.getStackTrace());
    }

    Element root = configFileName.getRootElement();
    if (null == root) {
      throw new ReadConfigException("Got no root element in config file");
    }
    else {
      this.cloudTemperature = Double.valueOf(root.getChild("temperature").getValue());
      this.cloudDensity = Double.valueOf(root.getChild("density").getValue());
      this.cloudGrainsNumber = Double.valueOf(root.getChild("number_sites").getValue()).intValue();
      this.cloudAccretionBatchSize = (Double.valueOf(root.getChild("accretion_step").getValue()))
          .intValue();
      this.iterationsLimit = (Double.valueOf(root.getChild("iteration_limit").getValue()))
          .intValue();
      this.cloudGrainSize = Double.valueOf(root.getChild("grain_size").getValue());

      this.cloudPopulation = new HashMap<String, Specie>();
      this.parsePopulation(root.getChild("population").getChild("gasphase").getChildren());

      this.reactions = new ArrayList<Reaction>();
      this.parseReactions(root.getChild("reactions").getChildren());
    }
    updateAccretionRates();
    initialize();
    setChanged();
    notifyObservers("model : config " + configName + " parsed at " + this.calendar.getTime());
    setChanged();
    notifyObservers(getCloudPhysicalParameters());
    setChanged();
    notifyObservers(getcloudReactions());
    setChanged();
    notifyObservers(getcloudPopulation());
    this.ready = true;
    writeSnapshot();
  }

  private Object getcloudPopulation() {
    CloudPopulation pop = new CloudPopulation();
    for (String name : this.cloudPopulation.keySet()) {
      pop.add(this.cloudPopulation.get(name));
    }
    return pop;
  }

  private CloudReactions getcloudReactions() {
    CloudReactions rStore = new CloudReactions();
    for (Reaction r : this.reactions) {
      rStore.add(r);
    }
    return rStore;
  }

  private CloudPhysicalParameters getCloudPhysicalParameters() {
    return new CloudPhysicalParameters(this.cloudTemperature, this.cloudDensity,
        this.cloudGrainsNumber, this.cloudGrainSize, this.iterationsLimit,
        this.cloudAccretionBatchSize);
  }

  public void run() {

    this.stopSimulation = false;

    boolean doLogging = false;
    if (null != this.fileWriter) {
      doLogging = true;
    }
    try {

      updateAccretionRates();

      makeUpExperiment1();

      for (int i = 0; i < this.iterationsLimit; i++) {

        StringBuffer sb1 = new StringBuffer(500);
        StringBuffer sb2 = new StringBuffer(500);

        sb1.append("## iteration : " + i + "\n");

        Long startTime = System.currentTimeMillis();

        TreeMap<String, Specie> acrretionSample = cycleAccretion();

        Long endTime = System.currentTimeMillis();
        String timeSpan = getTimeSpan(startTime, endTime);

        sb1.append("  model: accretion cycle finished in " + timeSpan + "\n");
        sb1.append("  accreted " + this.cloudAccretionBatchSize + "\n");
        for (String str : acrretionSample.keySet()) {
          sb1.append("  " + str + " " + acrretionSample.get(str).getAccumulatorValue() + " i.e "
              + Double.valueOf(acrretionSample.get(str).getAccumulatorValue())
              / this.cloudAccretionBatchSize + "% \n");
        }
        sb1.append("  ----\n");
        setChanged();
        notifyObservers(sb1.toString());
        if (doLogging) {
          System.out.print(sb1);
          this.fileWriter.write(sb1.toString());
        }

        startTime = System.currentTimeMillis();
        SurfaceCycle(acrretionSample);
        endTime = System.currentTimeMillis();
        timeSpan = getTimeSpan(startTime, endTime);
        sb2.append("  model: surface cycle finished in " + timeSpan + "\n");

        HashMap<String, Specie> grainPopulation = getGrainPopulation();
        sb2.append("  grain population: " + grainPopulation.size() + "\n");
        for (String str : grainPopulation.keySet()) {
          sb2.append("  " + str + " " + grainPopulation.get(str).getAccumulatorValue() + "\n");
        }
        setChanged();
        sb2.append("  ---- ---- ---- ----\n");
        notifyObservers(sb2.toString());
        if (doLogging) {
          System.out.print(sb2);
          this.fileWriter.write(sb2.toString());
          this.fileWriter.flush();
        }

        if (this.stopSimulation) {
          break;
        }
      }
    }
    catch (IOException e) {
      // TODO Auto-generated catch block
      e.printStackTrace();
    }
    tearDown();
  }

  private void makeUpExperiment1() {

    // check and set velocities for species
    for (Specie sp : this.cloudPopulation.values()) {
      if (sp.getName().equals("H")) {
        continue;
      }
      if (sp.getName().equals("D")) {
        continue;
      }
      if (sp.getName().equals("CO")) {
        continue;
      }
      sp.setAccretionProbability(0D);
    }
    Double normalizer = 0.0;
    for (Specie sp : this.cloudPopulation.values()) {
      normalizer = normalizer + sp.getAccretionProbability() * sp.getVelocity();
    }
    for (Specie sp : this.cloudPopulation.values()) {
      sp.setAccretionProbability((sp.getAccretionProbability() * sp.getVelocity()) / normalizer);
    }

    // for (Specie sp : this.cloudPopulation.values()) {
    // System.out.println("Set " + sp.getName() + " " + sp.getAccretionProbability());
    // }
  }

  public boolean modelReady() {
    return this.ready;
  }

  public void setPhysicalParameters(CloudPhysicalParameters pp) {
    this.cloudTemperature = pp.getTemperature();
    this.cloudDensity = pp.getDensity();
    this.cloudGrainsNumber = pp.getSitesNumber();
    this.cloudGrainSize = pp.grainSize();
    this.iterationsLimit = pp.getIterationsLimit();
    this.cloudAccretionBatchSize = pp.getAccretionStep();
  }

  public void setLogFile(File file) throws IOException {
    this.fileWriter = new BufferedWriter(new FileWriter(file));
    setChanged();
    notifyObservers("File " + file.getAbsolutePath() + "selected for logging.");
  }

  public void tearDown() {
    if (null != this.fileWriter) {
      try {
        this.fileWriter.close();
      }
      catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }
    }
  }

  public void update(Observable o, Object arg) {
    if ((o instanceof Controller) && (arg instanceof Boolean)) {
      this.stopSimulation = (Boolean) arg;
    }

  }
}
